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ABSTRACT 


Estimating the frequencies of signals in a noisy environment has numerous 
applications in digital signal processing. In December 1980, Golub and Van 
Loan proposed a spectral estimator called the Total Least Squares ( TLS ) 
technique which is a refinement of the Least Squares ( LS ) technique. In this 
thesis, we first describe the origin of the TLS technique and its applications to 
frequency estimation. Furthermore, we present a numerical implementation for 
resolving two damped / undamped closely-spaced sinusoidal signals in white 
noise. 

Next, we introduce TLS extensions such as the Constrained Total Least 
Squares ( CTLS ) technique and the Linear Constraint Total Least Squares 
( LCTLS ) technique. The CTLS addresses the case where the noise 
components are related and the LCTLS addresses the case where one desires 
to resolve between two narrowband signals close in frequency, one of which is 
known. Finally, we present a numerical implementation of the Recursive Total 
Least Squares ( RTLS ) technique and apply it to the case of a signal with a 


fixed frequency together with a signal with a time-varying frequency. 
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I. INTRODUCTION 


Estimating the frequencies of signals in a noisy environment has numerous 
applications in digital signal processing. Various techniques have been proposed 
to solve that problem such as the Pisarenko algorithm [1, p. 623] and the 
MUSIC algorithm [1, p. 626]. These two algorithms adopt the eigen- 
decomposition technique to obtain the solution from the partitioned signal and 
noise subspaces. In this thesis, we present a subspace technique for frequency 
estimation called the Total Least Squares algorithm. Two extensions, the 
Constrained Total Least Squares algorithm and the Recursive Total Least 
Squares algorithm, are also introduced. 

In Chapter II, we first introduce the ordinary Least Squares ( LS ) 
algorithm [1, p. 518]. The LS algorithm is used to solve the overdetermined set 
of linear prediction equations Ax=b, where A and b are mxn and mx1 matrices 
depending on data, and x is an array of parameters. In other words, the LS 
algorithm computes the solution xyg so that the Euclidean norm ll Ax-b ll is 
minimum. Note that the LS algorithm considers the errors or perturbations in 
the observation vector b only. However, in practice, we also need to consider 
the errors or perturbations occurring in the data matrix A. Thus, we consider 
the Total Least Squares ( TLS ) algorithm [2] which is an extension of the 
classical LS algorithm and has been applied to frequency estimation in signal 
processing to locate signals in the presence of additive noise. This algorithm 


can be viewed as a generalized LS algorithm and it has better performance 


than the ordinary LS algorithm. However, it is computationally more expensive. 








We obtain the TLS solution xq,s by simultaneously minimizing the 
perturbations AA and Ab which exist in the data matrix A and the observation 
vector b, respectively. . 

When there exists some linear dependence algebraic relationship among the 
correlated noise perturbation components in A and b, the TLS algorithm has to 
be modified to take this dependence into account. The modification leads to the 
Constrained Total Least Squares ( CTLS ) algorithm [6]. To be more specific, 
the CTLS algorithm results from an unconstrained minimization problem over a 
small set of variables. It uses the fact that the elements in perturbations AA 
and Ab are algebraically correlated. By employing a new matrix F with a 
smaller dimension and a white noise random vector w, it minimizes the 
Frobenius norm of AC to derive the solution Xcpz5. Next, we introduce the 
LCTLS algorithm [10]. This algorithm can be considered as an extension of 
the TLS algorithm and a special case of the CTLS algorithm. It is used to 
resolve closely-spaced frequencies with a given known frequency component. 
An example and a summary of the main steps needed to be implemented the 
LCTLS are presented. 

Chapter III presents a numerical implementation of the TLS algorithm for 
resolving two closely-spaced damped / undamped sinusoids. We investigate the 
effects of choosing different signal phase angles, signal-to-noise ratios ( SNR ) 
and prediction filter orders. In addition, we explore the estimation behavior of 
the algorithm due to incorrect singular subspace partitions. 

In Chapter IV, we introduce the Recursive Total Least Squares ( RTLS ) 
algorithm [11]. The algorithm is used in this thesis to estimate the frequencies 


of two complex sinusoidal signals in white noise. We adopt a spherical 








subspace updating [12] technique and use an eigendecomposition method to 
obtain the solution xpyys. Four parameters are varied: the signal-to-noise 
( SNR ) ratio, the prediction filter order, the forgetting factor and the range of 
the frequencies involved. We also study the algorithm performances obtained 
by using the signal or the noise subspace to estimate the frequencies. Last, 


Chapter V presents conclusions. 


Il. THE TOTAL LEAST SQUARES TECHNIQUE 
AND ITS EXTENSIONS 


A. INTRODUCTION 

Resolving two closely spaced sinusoidal signals in a noisy environment is a 
problem of importance in digital signal processing. This problem becomes 
rather difficult when the number of data samples is small or the signal-to-noise 
ratio ( SNR ) is low. Here, we adopt the Total Least Squares ( TLS ) approach 
to solve a Linear Prediction Equation ( LPE ) to improve the resolution. Note 
that the performance of the algorithm deteriorates when SNR is low, say below 
3dB. 

The TLS technique is the generalization of the ordinary Least Squares 
(LS ) technique [1] and it is often used to solve an overdetermined set of linear 
equations obtained from noisy observations. When there are more equations 
than unknowns, the system Ax=b is said to be overdetermined. The TLS 
method takes into account the errors that occurs in the data matrix A as well 
as in the observation vector b and try to attenuate these noise effects on both 


sides simultaneously. Further details are presented later in this section. 


B. LEAST SQUARES ( LS ) ALGORITHM 
1. Introduction 
First, recall the ordinary Least Squares ( LS ) technique. For the LS 
algorithm, the data matrix A is assumed to be free of error and all errors are 


confined to the observation vector b. However, this assumption is not always 


realistic since many errors ( perturbations / interference ) like sampling errors, 








modeling errors and instrument eure: etc. may result in inaccuracies or 
uncertainties in A. — | 
2. . Algorithm Derivation 
We obtain the LS solution x by minimizing the Euclidean norm 


| Ax - b Il. Thus, it is equivalent to the following expression: 


min||Ab|| | (2.1) 
x,Ab . 
subject to 
Ax=b+Ab (2.2) 
which leads to 
X= ( AMA) A®b = A*b (2,3) 


where A* is the pseudoinverse of A, A®A is the correlation matrix. The 
superscript notation 'H' denotes the complex conjugate transposition. Note 
that the LS algorithm works well when considering the noise only in the 


observation vector b. 


C. TOTAL LEAST SQUARES ( TLS ) ALGORITHM 
1. Introduction 
First, we represent the forward linear prediction equation Ax=b in a 


matrix form: 


Eo Er ots Epa If Xp Ep 


51 52 2 SP a - Spel (2.4) 


Enpa yep Eno dh oN-1 : 


iF is the observed noisy signal, P is the prediction filter order and is usually 
larger than the number M of exponentials but smaller than the number N of data 
samples. 

If the observed data is noisy, from a statistical viewpoint the LS 
solution will no longer be optimal because it can be biased and suffer 
covariance from the accumulation of noise errors in AMA. To alleviate this 
problem, we adopt the Total Least Squares ( TLS ) technique [2]. It has been 
shown to be equivalent to the minimum norm method [3] and it is used to 
diminish the bias accrued in AMA or to remove the noise components by using 
a perturbation analysis on A and b of the smallest 2-norm that makes the 


system equations consistent. The TLS solution can be obtained by minimizing 


WAC I = i [AAlAb] I (2.5) 
subject to 
(A+AA)x=b+Ab (2.6) 
which leads to 
Xorg = ( AMA - O1)'AMD | (2.7) 


where AA and Ab are two independent noise perturbations and o” is the 


minimum eigenvalue of [ Alb }2[ Alb }. 











Note that we use a zero mean white identically iidependint 
distributed ( i.i.d. ) Gaussian noise in our analysis of the TLS algorithm. When 
the noise is not white but colored and still uncorrelated with the signals, a 
whitening filter ought to be applied. This kind of filter is capable of transforming 
a stationary discrete-time non-white process at the input into a white 
uncorrelated data sequence at the output. The whitening procedure takes the 
non-white noise vector w and computes R=E[ww*]=ZZ"*, using a Cholesky 


decomposition. This leads to the definition of a new white noise vector 5 where 
6 = Zw. (2.8) 


The whitening procedure can be performed by using a prediction-error filter 
which offers a linear prediction equation same as Ax=b. Basically, the 
prediction is dependent on the presence of correlation between adjacent 
samples of the input non-white stochastic process. Then, we can successively 
decrease the correlations by increasing the prediction-error filter order high 
enough until it ultimately gets to a point at which the output sequence is 
composed of ( white ) uncorrelated samples [4, p. 216]. 
The TLS solution X71, is obtained by using information obtained from 
the singular value decomposition ( SVD ) of the matrix C=[ Alb ] with a 
dimension (N-P)x(P+1), i.e., 
C = Uxv" (2.9) 


where U and V are unitary square matrices with a dimension (N-P)x(N-P) and 
(P+1)x(P+1), respectively; Z is a rectangular diagonal matrix with a dimension 


(N-P)x(P+1) consisting of the singular values in decreasing order with the largest 





one on the upper-left-handed corner. The matrices U, V and Z can be 


partitioned as 


U=[u,lug!---luygp] | (2.10) 
Ve=lv,ivyt---l vp.) ] (2.11) 
and X= diag (0),--+, Op Op,; )- (2.12) 


The columns of the matrices U and V are the unitary eigenvectors of AA# 
and AMA respectively. 
2. Algorithm Derivation 
The TLS algorithm described above can be summarized into the 


following equations [5] : 


({ Alb ]+[AAIAb J) (*) =0 2.13) 
or 
oakeist 7 (2.14) 
where 
C=[Alb ], AC =[ AAIAb ], z= fe (2.15) 


From the homogeneous equation (2.13), we see that the TLS problem 


finds the minimum norm of the perturbation matrix AC such that ( C+AC ) is 





rank deficient. In other words, the TLS solution can be formulated as seeking a 


solution vector x [6] to the following problem 


infACIZ 
min| Ir. 


(2.16) 
subject to 
x 
(c+ac)(")} =0 | (2.17) 
where 
_ ACE = tr{AC* AC}, | (2.18) 
ll - Ip denotes the Frobenius norm and" tr" is the trace of a matrix. 
Next, the matrix C is decomposed in the following form 
| x, O] Ve 
~ 2. 
com off SRE] 


where Uj, 21 and V, correspond to the signal subspace; U2, Z2 and V2 

correspond to the noise subspace. We separate the signal subspace from the 
noise subspace by using the number of sinusoids. Two different cases need to 
be considered. The first one considers the case where the singular values are 
all distinct. The second one considers the case where some of the singular 
values are multiple. For the first case, we usually assume that the last element 
Op,; in the diagonal matrix B is the smallest singular value and it corresponds to 


the last vector Vp,; in the unitary matrix V. Then we get 





1 V(P +1), 
XTLs ="-77 : (2.20) 


(vps Ve V(Pat), 


where (vp,;)p,, is the last element of the last singular vector vp,; and is 
assumed to be non-zero for normalization. If (vp,)p,, is equal to zero, a TLS 
solution does not exist. For the second case, we assume that the P-M+1 


smallest singular values are equal, i.e., 


Oy 202 2-2 Oy > OMe = °° Spy, | (2.21) 


and set V2 to be a column partition of V where 
V2=(Vyap tts Vp). (2.22) 


or 


cl 
V>= & . (2.23) 
2 


where cT represents the first row of the matrix V>. 
Next, we introduce the Householder ( reflection ) matrix Q which is 
used to solve the TLS problem. It has the property to zero out specific entries 


in a matrix or vector. Usually, it is defined as [4, p. 428]: 


* 


VV 
* 
vv 





Q=I,-2 (2.24) 
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where 
I is an mxm identity matrix, 
v is an mx! complex vector which corresponds to the complex vector v; in 
V> defined in (2.22) and has the Euclidean norm 


i 
Ivl=(v"v)? (2.25) 


The vector v can be defined below for analysis: 





vec -aep = Vi=cl-atef (2.26) 

on =t2Icl = a=t-Pfe| (2.27) 
|cp| cp 

eb=(0, ---, 0, Ip => Vp =cp-& (2.28) 


For any given vector, we can construct a Householder matrix so that all the 
energy is compressed into a selected partitioned column vector [7]. This leads 


to the following TLS solution: 


0 « 2.29) 
¥:0-| > 


c'Q =[0,---,0,a]. . (2.30) 


Now we express VQ to extract the last column, which is used in the 


solution: 


ee EE.” 





and let 
a 
Vod5 -[* . (2.32) 


Note that 


viy =(c7-aTep \(c’-a'ep) 





= lel? - a*cp - acp +a? G a cp = acp ) 
* 
c 
= lel? -20°cp +|}af? oO =t—* el 
lep| 
2 leer lepl 2 A, ene 
=|cl" - lle|+.— sell ( “lal? = ao ) 
Teel P| |cp| 
cpe ; 
= {ier : otet| (s lla? = lle? = cHe) 
lep| 
= 2(ccr” -a"cp). ; (2.33) 


The last column of the matrix Q is given by 





"at a) 


ae Cp 
P 20° 2a" (ate) Cp ) 


Pen) ies -0"ep ) 
* * 
Cpe” -Cpa ep - ac +O0Q% Cp 
° a" (a- cp) 


c’ (a-cp) 
(a-cp) @. (2.34) 
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Thus, 


' Vie 
y= Voq >= ae nae | (2.35) 
a 
and 
y Vc" Vic 
Xq co =e S-—e =- 2.36 
TS a Oo cle eS?) 


Equation (2.36) corresponds to the noise subspace case since V; and ¢ are 
obtained from V2 which is described in (2.19). Similarly, we can derive the 


solution vector from the signal subspace. First, we use the identity 


H .H H 
Bee Maley (2.37) 
V; V, ¢c V, . 
which leads to 
gig+cHe=1 => cle=1-gg 
Vig + Ve =0 => Vic=-Veg . (2.38) 


Next, we substitute (2.38) into (2.36) and obtain the alternative TLS solution 
corresponding to the signal subspace case: 
Viz" 


= 2.39 
ca (2.39) 


XTLS 


where V; and g are obtained from V, which is also described in (2.19). 


i3< 





D. CONSTRAINED TOTAL LEAST SQUARES (CTLS) 
ALGORITHM 
1. Introduction . 
When the noise components contained in the augmented matrix C 
are algebraically correlated or linearly dependent ( so called constraints ), the 
TLS technique is no longer an optimal frequency spectral estimator. In this 
case, we have to reformulate the problem to reduce the dimensionality of the 
disturbances. Thus, the TLS needs to be reformulated to reduce the 
dimensionality of the noise entries. Thus, the Constrained Total Least Squares 
( CTLS ) algorithm, proposed by Abatzoglou et al [6], can be viewed as the 
natural extension of the TLS algorithm. In addition, the CTLS can be shown to 
be equivalent to the Maximum Likelihood ( ML ) algorithm. The CTLS solution - 
can be obtained from a constrained state-space-parameter ML estimator 
together with a Newton algorithm [8]. | 
2. Algorithm Derivation 
In order to remove the correlated noise perturbation components 
among A and b, we derive the CTLS algorithm from an unconstrained 
minimization problem over a small set of variables. Let us redefine AA and Ab 


both as linearly constrained. 


AA =[ F,w!Fow!--:|Fpw] | (2.40) 


Ab = Fpi.W (2.41) 
Then the augmented error matrix AC is formed by 


AC =[AAIAb]=[F,wIF,wl---!Fpw|Fpyw] (2.42) 
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Or 
AC =[AC,!---|ACp| ACp,; ] = Fw, (2.43) 


F=[FyIF,l-s-IFplFeal (2.44) 


where F is a matrix whose elements are related to those of independent 
variables AC and w is a zero mean white noise random vector of minimal 
possible dimensionality. The reason for w holding the minimum dimension is to 
increase the number of algebraic relationships the entries of AC can satisfy and 


to accomplish the minimization faster. So the CTLS solution can be formulated 


as follows: 
min||w|? (2.45) 
w,Xx 
subject to 
r . 
(C+AC ) ( }- 0. | (2.46) 


Next, we consider the quadratic constraint equation in (2.46) which is expanded 


c{*]+ac{* l= 0 ie 2.47 
(i}+ac(i}- ean 


into: 
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where 


x | x 
a eo] 


Pp P 
= > x,Kw-F4w {y xX - Foy \ me, (2.48) 
i=l j=l é 
Let 
= 
H, = ¥xF = Fou; ’ : (2.49) 
i=l 
then we have 
. +c{ * 2.50 
Cc e +H,w=0 = w=-H,C 4 (2.50) 


where H7¥ is the pseudoinverse of Hx defined as 
Hi = Hi(H,H;) | (2.51) 
For any w, the solution vector x can satisfy (2.45) and (2.50); i.e., we have 
min|w/? 2 (*) eiitol), (2.52) 


Next, let us define the analytic function J(x) 


* 


J@e= (*) c'(az) #ic{*) (2.53) 
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X=(Xpe-+, Xp). 
J(x) can be simplified to the following expressions by substituting (2.51) into 


(2.53). | 
xX * a \-L ae *\ % 
sx)=(7) Cc H, (HH; ) 'H; (HH) ce 


2 (*) c (HA ) HH, |(H,H, Ng c{”) 
= Ge c'(H,H,)'¢{ *} | (2.54) 


Hence, we conclude that the CTLS solution is obtained by minimizing 
J(x). From (2.45) and (2.52), we try to find out the minimum of J(x) with a 
respect to x. It is necessary to compute the first-order complex partial 


derivative or complex gradient and set the derivative to be zero; namely, 


ay _( 23 ot). % 
21 (38 3h) 0, a) 


Due to the complex data, Abatzoglou et al [6] defined a complex 
version of the Newton method which uses the second-order complex partial 
derivative of J(x) to update the values of x. The Newton iteration ( recursion ) 


is briefly described by 


x(n-+1)=x(n)+(ABTA -B) (a-AB“a) (2.56) 
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: | hg 
= oy = oot) _ complex gradient of J 
Ox \0x, OXp 
2 2,T 
= (3 a = unconjugated complex Hessian of J 


o7J 
B = = + . 
Axx conjugated complex Hessian of J. (2.57) 


Hence, by referring to the derivation in [6], we obtain the following 


closed-form solutions for a, A and B: 


a=(h*B)" 
A=-—F" H¢(H,H*)" B-(F’ H¢(H,H’)' B)" 
B=[B’(H,H,)' B]’ + F [H, (H,H,)” H, -IJF (2.58) 
where 
* -] x 
h=(H,H,)C 1 
B=CIp.ip —[F,Hjh|-- |FpH,h] 
G =[F,h|---|Fph]. : (2,59) 
The Newton method is implemented by starting with an /terative 
Quadratic Maximum Likelihood ( IQML ) algorithm [9]. This algorithm solves 
a quadratic constraint minimization problem and derives the solution vector 
from the coefficients of the linear prediction polynomial which correspond to the 


parameters of the signals in a fast convergent way. We summarize its main 


steps below: 
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(a).initialization : set x(0)=xo=1 at n=0, x is viewed as the vector composed of 


the prediction polynomial coefficients and obtained using Newton iteration 


T 
x =(Xp.X15"7"Xp] 


(b).compute : 
N 
Zy => YP (XX); 
jel 
where . You Yo 
Y hie i Yi 
Ym-1 Ym-p-2 Ym-p-1 


c).solve the Quadratic Minimization Problem : 


‘ H (n) 
min {x Zy’X 
‘ main ( (nt+1)“Y “(n+1) ) 


(d).set n=n+1 for real-time iteration 


(e).check the convergence: II x(n-1) - x(n) Il<e 


(2.60) 


(2.61) 


-(2.62) 


(2.63) 


(2.64) 


(f).identify the roots of the prediction vector polynomial x which represent the 


parameters of the received signal. 
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E. LINEARLY CONSTRAINT TOTAL LEAST SQUARES 
(LCTLS) ALGORITHM 
1. Introduction 
In this section we introduce a special case of the CTLS algorithm 
called the Linearly Constrained Total Least Squares ( LCTLS ) technique. The 
LCTLS technique is used to locate and resolve closely spaced frequencies in 
the presence of a given known frequency component [10]. 


Recall that the TLS method solves the unconstrained problem : 
(A +AA)x=b+ Ab. (2.65) 


Using (2.13) to (2.18), we compute the minimum norm of the total perturbation 
matrix AC which makes ( C + AC ) rank deficient and generates a vector z 


defined from a null space. Equation (2.65) has a solution provided 
(b+Ab)€ Range (A +AA ). (2.66) 


In this case, 2,4;#0 and we can solve for x. If zn1:=0, the TLS problem has no 
solution, In addition, the non-zero last element z,., is used to normalize the 


vector Z; 1.e., we have 


{2 |.[ * a oe a : 2.67 
(Z(G) = (3) ca) 


2. Algorithm Derivation 





To solve the LCTLS problem, we first look at the constraint equation: 


ypx=T (2.68) 
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where , x and I’ have dimensions px(n+1), (n+1)x1 and px1, respectively. 


Equation (2.69) can also be expressed as below: 


(we ri(7}=19 []z=0. (2.69) 


Let us introduce a physical example to define the matrix [ ‘¥ I]. First, we 
consider a signal €(n) composed of two complex sinusoidal signals in white 


noise: 


C(n) = Aye" + A,e™™ + w(n), n=1,--,N (2.70) 


where C(n) is the observed noisy data; w(n) is the white Gaussian noise; A; is 
the amplitude; «; is the angular frequency and N is the data length. 
Assume that @, is known and 2 is unknown. Create a signal 


source vector A(w) defined as: 


AA(w) =[1, 3, e?, cP? 1. eM] (2.71) 


We make A(m) satisfy the constraint equation in (2.69), which leads to 


aH (2) =0. (2.72) 


The constraint forces e to be one of the roots of the prediction-error 


polynomial e(@) defined as: 


e(o) = m*o( ori (9,73) 


2h 





Then, the angle of the next closest root on the unit circle is chosen as the 


estimate of the unknown angular frequency @». 


Next, considering the TLS equation given in (2.65) and constraint 


2 OF ay | 
A+AA b+Ab (7,)- eR _. G74) 


Next, using a SVD partition method to obtain the solution. To ensure the 


given in (2.69) leads to 


constraint equation in (2.75) to be satisfied, we have to constrain the vector z to 
fall in the null space of { ‘¥ I] which can be further factored out using OR 


decomposition into 


[W P]=[Q,0j]A7=0,A7 (2.75) 
H ; | 
Ate (* | (2.76) 
2 


where Q, is an upper triangular matrix with a dimension (pxp) and AX is a 
unitary matrix with a dimension (n+1)x(n+1). By comparing the matrix 
dimensions from (2.69) and (2.75), we see that z should be the product of the 
orthonormal basis A, with a dimension (n+1)x(n+1-p) for the null space of 


{ Y I] and an arbitrary non-zero vector 7 with a dimension (n+1-p)x1, ie., 


Z= AN . (2.77) 
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Next, define the matrix C as 

C=CA,Aq | | (2.78) 
with 

C =[AIb]. | (2.79) 


The SVD of the matrix C is performed and C is partitioned into its signal and 


noise subspaces, which leads to: 


~ «~ «0H n+l-p ~ «~H 


C=UEV = ¥ 9, u,v; (2.80) 
i=l 
or 
zy, 0 OlvE 
é-[¥, U, Us| 0 x, Of VE 
0 0 | v8 © 2.81) 
with 


V1 :1-dimensional unconstrained part of the signal subspace 
V2_ : s-dimensional noise subspace 

V3 : p-dimensional constrained part 

Dy cineoiees GO) 

2, = Snap I,, S=n+l-p-r 


23 =0. 
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Using the fact that 


A, AU +A, AE =I, | (2.82) 
we get 


C=CAAS+CAAS (2.83) 


Dowling et al [10] have suggested to enforce sphericity in the 
concerned partitioned subspace to improve the frequency spectral estimation 
performance. The term " sphericity “ means that the singular values in the 
signal and noise spaces are replaced with their mean values. For example, 
‘assume that there exists ( s=n+1-p-r ) equal singular values in the noise 
subspace V2: ie., the noise singular values are replaced with the variance of 
the random noise. By partitioning the truncated unitary matrix V below and 
doing the same operations as in (2.37) and (2.38) 

v-|¥, V, ¥,]-|" V2 | (2.84) 


gt ct ft 


we obtain the minimum norm solution corresponding to the noise subspace: 








Further using the following identity : 


H 
V; 
gi cH fH ae 
a cake =I, (2.86) 
1 2 3 f V3 . 


which leads to the following two equations: . 
gloscHe+ftf=1 = cHe=1-gig-fMFf 


VigtV,c+V,f=0 = V,c=-Vig-Vaf (2.87) 


we obtain the minimum norm solution corresponding to the signal subspace: 


Vie +V;f- 


XLCTLS = 1-oWe- fhe oe - fhe (2.88) 


Note that the minimum LCTLS solution should be based on the subspace V, or 


ki V; | whichever one has the smaller dimension. 
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WI. TOTAL LEAST SQUARES TECHNIQUE 
IMPLEMENTATION 


A. INTRODUCTION 

First, we study the performance of the TLS algorithm to estimate the 
frequency of sinusoidal signals in noise. We assume the received signal x(n) 
consists of M uniformly spaced signal data samples corrupted by a complex 
white Gaussian noise w(n). The received signal is expressed by the following 
equation and is observed over a duration of N seconds, assuming the sampling 


frequency f,=1Hz : 


x(n)= 5 a," 4.w(n), n=0, 1-,N-1 | (3.1) 
. kel 
"where 

{ a, } : aset of signal amplitudes 

{ f, } : a set of true signal frequencies 

{ x(n) } : measured complex-valued data samples 

{ w(n) } : random complex-valued noise samples and uncorrelated 

with signals 

Real and imaginary part of w(n) are assumed to have a variance 07. 

We are interested in estimating the signal frequency parameters f,; the 
parameters a, are chosen to have unit magnitude for the sake of simplicity. 
Note that we use the singular value decomposition ( SVD ) technique to solve 
the linear prediction equation Ax=b which matrix form is given in (2.4). The 


TLS implementation is expected to decrease the noise effects occurring both in 
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the data matrix A and in the observation vector b simultaneously. First we 
choose the order P of the prediction-error filter and solve a P-degree prediction- 
error filter polynomial. Next, we determine the signal parameters from the M 
zeros which are the closest to the unit circle. The other P-M zeros are called 
extraneous zeros; they are approximately uniformly distributed inside the unit 
circle. Note that the choice of optimal prediction-error filter order is important 
since if the filter order is chosen too low, the closely-spaced sinusoids can not 
be resolved. In addition, if the filter order is chosen too high, the noise zeros 
may be mistaken for the true signal zeros as the noise zeros become very close 
to the unit circle. Thus, we have to select the " optimal " prediction filter order 


carefully. 


B. ALGORITHM IMPLEMENTATION 


? ; . 119 
The parameters used in the simulation are =F o,=0; f; and f, represent 


the frequency locations; w(n) are independent identically distributed ( i.i.d. ) 
white complex Gaussian noise samples with zero mean and variance 0? for 
both real and imaginary parts. The signal-to-noise ratio ( SNR ) of the received 


noisy signal is given by the peak amplitude of sinusoids to noise level and 


defined as 10slog( 57 MB. The sampling frequency is set to be 1Hz. The 


frequency separation is chosen smaller than the reciprocal of the observation 
time in order to resolve these two close-together frequencies and raise the 
resolution. Simulations are performed using an ensemble of 30 independent 
trials. For each trial, a data block of 50 ( N=50 ) data samples is used to 


resolve two ( M=2 ) sinusoids. The frequency estimates are computed by 


pe 





choosing the two roots of the prediction-error filter transfer function or the two 
zeros of the prediction-error filter polynomial which are the closest to the unit 
circle in the complex z-plane. Different values of the filter order P in the range 


of 2~49 can be used to observe the variations of estimation. 


C. EXPERIMENTAL OBSERVATIONS 
1. Performance of the Algorithm 
We investigate the performance of the algorithm with a statistical 


quantity called the effective standard deviation ( ESD ) defined below : 


1 Tyna 72 | 
ESD = — > lf,-f (3.2) 
Tin 








where | 
T : the total number of independent trials 
f : the estimated frequency at the it realization 
f : the true signal frequency. 
The ESD gives a measure of the distance between the true and the estimated 


frequency locations. In this simulation, we first use f,-0.23Hz while f, =0.27Hz. 


The data including =F, $,=0, P=8, N=50 and 30 independent trials are used 


to compute the ESD. The selection of the two signal zeros is made by creating 
two sector regions designated from the origin to point z,=0.3+j*0.7 and to point 
z= -0.3+j*0.7, respectively. Varying the SNR from 20dB to OdB, we observe 
that the sector region for f;=0.23Hz at SNR=0dB needs to be expanded by 
changing z,=0.4+j*0.7 and the sector region for f-=0.27Hz at SNR=1dB needs 
to be expanded by changing z,=-0.4+j*0.8. Here we display six different SNR 
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cases which are 20dB, 12dB, 5dB, 3dB, 1dB and 0dB. We also observe that the 
signal and noise zeros resolution becomes worse as the SNR becomes smaller, 
as shown in Figures 3.1.1~3.1.6. Figure 3.1.6 shows the worst resolution as a 
couple of noise zeros are on the unit circle and could potentially be considered 
as the signal zeros. In addition, Table 1 shows that as the SNR becomes 
smaller, the estimated deviations for both signal frequencies are getting larger. 
The above simulation results compare with those obtained in [5]. 
2. Effects of Signal Phase Changes 


Two sets of phases are chosen to investigate the effects of phase 
. T T qT ; 
changes: (1) oiea: $2=0, (2) 1=7 ar The filter order is P=8; the 


number of data samples points is N=50; 30 independent trials are used; the 
signal frequencies are chosen as f;=0.25Hz, f2=0.27Hz; the SNR values — 
considered are 3dB and 12dB. Figures 3.2.1~3.2.4 show that phase changes do 
not affect the resolution. 
3. Effects of Filter Order Changes 

Next, we consider the effect of changes in the prediction-error filter 
order. Recall that it is important to choose the smallest value of P which 
achieves a good frequency estimation performance. For this experiment, we 
choose f;=0.25Hz, fo=0.27Hz, SNR=12dB, N=50, and 10 independent trials are 
used. Figures 3.3.1 to 3.3.6 show the effects of changing the filter order. When 
P=2, we have very poor resolution, as shown in Figure 3.3.1. Performance 
slightly improves when P=4, as shown in Figure 3.3.2. Note that increasing the 
filter order P improves the resolution of the estimated frequencies. For 


example, when choosing P=8 and P=18, we observe that two sinusoidal signals 
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are located exactly at their desired positions on the unit circle and the 
extraneous zeros are uniformly distributed inside the unit circle, as shown in 
Figures 3.3.3 & 3.3.4. However, if P is chosen too large, P=28 or P=38 for 
example, we observe some of the P-M extraneous zeros of the prediction 
polynomial are getting very close to the unit circle which make the distinction 
between signal and noise zeros difficult to make, as shown in Figures 3.3.5 & 
3.3.6. 


4. Effects of Errors in the Subspace Partition Estimate 


It is possible to express the TLS solution in the sense of minimum 
norm. By applying the SVD technique to the linear prediction equation, we can 
partition the unitary matrix V ( for handling complex data ) into two subspaces: 
the signal subspace and the noise subspace. If the partition is performed 
correctly, we can obtain a good zeros resolution by using a moderate filter 
order P and high enough SNR, choosing P=8 and SNR=12dB for example, as 
shown in Figure 3.3.3. However, if the partition is not performed correctly; i.e., 
the numbers of column vector of the signal subspace and the noise subspace 
are not estimated properly, we observe that the extraneous zeros of the 
prediction polynomial are randomly distributed. The more the signal and the 
noise subspaces are incorrectly chosen, i.e., the worse the partition is, the 
worse the resolution to identify the true sinusoidal signal locations becomes, as 
shown in Figures 3.4.1 & 3.4.2. Figure 3.4.2 displays the worst degradation of 
the performance since these two subspaces are mixed up by 6 column vectors 


under the condition of SNR=12dB, P=8, N=50 and 10 independent trials. 
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5. Application to Damped Signals 
In the following, we consider the case of a sum of two ( M=2 ) 


damped exponential sinusoids described as follows : 


x(n)= ef Oumrslenhn+9,)) 4 g(-anri2nn+o,)} y(n), (3.3) 


n=0,1,---,N-l 


where angular frequencies are f,=0.25Hz, f,=0.27Hz; the signal phases are 
fixed at eT, ,=0; the two damping factors are chosen @)=0.8, a@,=0.9; the 


number of data sample is N=50; the filter order is P=8; 30 independent trials are 
used for the simulation, and w(n) is still defined as the complex white Gaussian 
noise. We observe that it is quite difficult to identify the exact locations of these 
two damped sinusoidal signals at low SNR, for example 5dB, or at high SNR, 
for example 20dB, as shown in Figures 3.5.1 & 3.5.2. Even if we increase the 
filter order up to P=18, the resolution has not been improved still, as shown in 


Figures 3.5.3 & 3.5.4. 
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IV. RECURSIVE TLS TECHNIQUE 
IMPLEMENTATION 


A. THEORETICAL BACKGROUND 
1. Introduction 

As discussed before, the Total Least Squares ( TLS ) algorithm has 
been shown to be more accurate than the Least Squares ( LS ) algorithm when 
the errors or noise perturbations exist on both sides of a linear prediction 
equation Ax=b. We now introduce an adaptive algorithm called the Recursive 
Total Least Squares ( RTLS ) algorithm. In this algorithm, we perform the 
updates by tracking the smaller one of the signal subspace or the noise 
subspace partitioned from a unitary matrix containing complex data. It is 
shown that the RTLS algorithm works better than the Recursive Least Squares 
( RLS ) algorithm [4, p. 480] and the Least Mean Squares ( LMS ) algorithm 
[4, p. 302] in terms of tracking property as well as steady-state tap-weight error 
norms [11]. Furthermore, the updates can also be done using a computationally 
efficient non-iterative subspace updating technique [12]. 

2. Algorithm Derivation 

First, let us give a brief review of an adaptive process. An adaptive 
process involves the automatical adjustment of a set of tap weights of a 
prediction-error filter; i.e., it obtains an estimate of a desired response which 
results from the inner product of a set of tap input and its corresponding set of 
tap weights produced by the adaptive filter. Then, it compares this obtained 


estimate with the actual values of the desired response and generates an 
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estimation error. The estimation error is fed back to the adaptive process and 
stimulates the prediction filter in a recursive way. Therefore, this is a closed- 
loop operation and the final goal of the adaptive process is to make the output 
error reach its minimum value. There are two good examples: the LMS 
algorithm and the RLS algorithm. 

In the TLS algorithm, we apply the singular value decomposition ( SVD ) 
technique to the augmented matrix spanned by the data matrix A and the 
observation vector b. However, in the RTLS algorithm, we apply an eigen- 
decomposition technique to a specific augmented matrix C(n) defined below to 
get the desired subspaces and determine the effective rank of C(n) from its 
singular values. In addition, we adopt a parameter @ called the exponential 


weighting factor, which is similar to ) ( forgetting factor ) in the RLS algorithm. 


Cay= DM [Alb] (4.1) 


where 


D® = ./diag(a”,...,0°). | (4.2) 


It is worthy to investigate the impact of the forgetting weight factor & on 
the tracking capability. The factor a is used to ensure that the data samples in 
the distant past are " forgotten " so as to afford the possibility of following the 
statistical variations of the observable data when the filter operates in a non- 
stationary environment. It is designated to be a positive constant close to, but 
smaller than one. Further speaking, the reciprocal of 1-a can be considered as 


a measure of " memory ". When @ is equal to 1, which corresponds to infinite 
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memory or most unforgettable case, the algorithm keeps all data samples and 
the estimated frequencies have low variance. On the other hand, when @ is 
approaching to 0, which represents the most forgetful case, the algorithm keeps 
only a few data samples and the estimated frequencies have high variance. 

As mentioned above, the RTLS algorithm applies the eigen-decomposition 
technique to the newly augmented matrix Cin). The estimated time-varying 


correlation matrix Rn) is expressed as 


Reg) = AR yay + (1- ©) cay ¥Gn) (4.3) 


where R(n) is defined as 


H 
Ra) _ Cay Cay, ‘ . (4.4) 


In this tracking ( updating ) method, we use the rank-one eigenstructure 
update recursively. In order to diminish the computational load, we use the 
sphericity property proposed by DeGroat [12]. Spherical updating forces both 
of the eigenvalues in the signal subspace and the noise subspace to be replaced 
with their mean values. Thus, we form two different spherical subspaces 
without changing the dimension of the diagonal matrix D of Ram). The spherical 
subspace updating might be viewed as the nucleus of the RTLS adaptive filter 
algorithm. We implement the eigendecomposition of the one-step past-time 
correlation matrix R(n-1) to get the updated decomposition at time n. Thus, we 


have the following updating 


H 
Rey = D-Day U a-1 (4.5) 
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: 
Ray = Vain | (4.6) 


where for all present time n 
H H = 
Uy Ueay = Ula Ui = 1 (4.7) 


and Dq) is still a diagonal matrix consisting of the eigenvalues. 

This updating algorithm is developed to track the smaller averaged 
eigenspace of the exponentially weighted correlation matrix Ra) under the 
condition of forced sphericity of the eigenstructure. To determine the right 
singular subspace of C(n), we need to obtain the eigenspaces of R(n) 
corresponding to the right singular subspaces of Cn). For later simulation, we 
get the vector y(n) extracted from the last row of matrix [ Alb ]. Next, for the - 
purpose of simplifying the derivation of the subspace updating, we neglect the 
factor ( 1-a ) and adopt a new time-varying correlation matrix R for the rank- 


one updating recursion 
Ro) = OR q) + Yan), (4.8) 


By substituting (4.6) into (4.8) and letting 
y=Ury : (4.9) 
we further get 


y = Ux (4.10) 
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and derive the following form 
R=U(aD+xx")UF (411 


Since both of the RTLS algorithm and the TLS algorithm use the SVD- 
based partition method, they have the same solution for the signal subspace and 
the noise subspace, respectively. Namely, we can rewrite (2.36) and (2.39) 
from Chapter II below: 
for the signal subspace, the solution vector is 
1-g%g : (4,12) 


’ 


V, 
XRTLS = 8 


for the noise subspace, the solution vector is 


Vic" 
XpTLs =— He (4.13) 





In simulations considered later in this chapter, we use two signal 
frequencies: one frequency is fixed and the other one is linear time-varying, so 
the smaller subspace is the signal subspace. We assume that these two 
complex sinusoidal signals have frequency components f,=0.27Hz and f; 


moving around 0.23Hz with a moving weight of 3% and 10%; the phase 
components are aT, ,=0; 10 independent trials and 10 recursive realizations 


are used in this simulation. Note that we have four parameters to vary: the 


filter order P, the exponential weighting factor a, the signal-to-noise ratio 


( SNR ) and the moving weight of the time-varying signal frequency. 
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3. Simulations 
First, we apply the RTLS algorithm and use the spherical subspace 
updating technique. Next, we compare the performance of the RTLS algorithm 
with spherical updating with that of the RTLS algorithm without spherical 
updating. A data length of 100 points and a window length of 50 points are 
used in the following simulations. In addition, we assign the forgetting factor 
to be one of the two values: 0.25 and 0.95. Next, we adopt one of the two 
moving frequency weights: 3% and 10%. We would also like to see when the 
resolution of signal zeros breakdowns by selecting various filter orders. Finally, 
by fixing «=0.95 and the moving weight=10%, we decrease SNR to 0dB and 
-5dB so as to observe the influence of low SNR on the resolution. The cases 
studied are listed as follows : 
a, Case (1) 
Fix SNR=12dB, P=8, moving weight=3% 
Change a=0.25, 0.95. 
b. Case (2) 
Fix SNR=12dB, P=8, imovine weight=10% 
Use o=0.95 
c. Case (3) 
Fix SNR=12dB, a=0.95, moving weight=10% 
Choose P=3, 5. 
d. Case (A) 
Fix P=8, a=0.95, moving weight=10% 
Decrease SNR=0dB, -5dB. 


a1 





B. EXPERIMENTAL OBSERVATIONS 
1. Estimated Signal Subspace Obtained with Averaging the 

Eigenvalues 

The signal subspace information is used to estimate the two signal 
frequencies. By using the data in case (1), for a=0.25, we observe that the 
locations of the zeros are randomly distributed. It is very difficult to distinguish 
signal zeros from noise zeros, as shown in Figure 4.1.1. The averaged 
eigenvalues of the spherical updated matrix are presented in Figure 4.1.2. If a 
is chosen much higher, say a=0.95, we observe that the resolution becomes 
very good and we can discern the fixed signal and the moving signal. It is 
noted that the noise zeros just congregate forming a small cluster and are 
distributed inside the unit circle, as shown in Figures 4.1.3 and 4.1.4. Hence, the 
experiments show that the case of «=0.95 has much better resolution since it 
keeps more data samples as shown in Figure 4.2.1. Next, we observe that it is 
quite obvious to distinguish the fixed frequency ( a slight shift ) and the time- 
varying frequency ( countable small circles ) using the data given in case (2). 
The averaged eigenvalues are also displayed in Figure 4.2.2. Thus, results 
show that the case with larger moving weight and higher forgetting factor 
( approaching to one ) has the best zeros resolution and tracking capability. 
Further speaking, the forgetting factor plays a more important role than the 
moving weight. By using the data in case (3), we observe that in the case of 
p=5 the " fixed " signal estimate shows some some variations. However, we 
are still able to differentiate signal zeros from noise zeros, as shown in Figure 
4.3.1. The performance gets worse when we use P=3, as shown in Figure 


4.3.2. Usually the filter order number ( P ) is chosen larger than the true signal 
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: : ae ae 
number ( M ), so we do not consider the case P=M here. Using the data in case 


(4), we observe that the resolution becomes worse at very low SNR since 
some noise zeros are very close to or on the unit circle and can be mistaken for 
true signal zeros, as shown in Figures 4.4.1 & 4.4.2. : 
2. Estimated Signal Subspace Obtained without "Averaging = 

Eigenvalues 

Next, we consider the estimated correlation matrix R(n) obtained 
without using averaged eigenvalues. Using the data in case (2), we obtain very 
good resolution of the fixed frequency and the time-varying frequency, as 
shown in Figure 4.5.1. The non-averaged eigenvalues of diagonal elements of 
the correlation matrix are displayed in Figure 4.5.2. Next, using the data in case 
(4) shows that the resolution becomes worse with very lower SNR, as shown 
in Figure 4.6 to be compared with Figure 4.4. Besides, from the plot of non- 
| averaged eigenvalues, we see the signal eigenvalues are still much higher than 
the noise eigenvalues as shown in Figures 4.6.2 & 4.6.4. 

3, Estimated Noise Subspace Obtained with Averaging the 

Eigenvalues 

Here, we use the spherical subspace updating method and choose the 
noise subspace to estimate the frequencies. The data are generated using the 
characteristics of case (2) given in section 3. The results are quite different 
from those obtained when choosing the signal subspace. Some of the noise 
zeros appear on the unit circle and can potentially be mistaken for signal zeros, 
as shown in Figure 4.7 to be compared with Figure 4.2. Next, using the data 
defined in case (3) earlier, we still obtain poor resolution, as shown in Figure 4.8 


to be compared with Figure 4.3. Next, we investigate the effects of the SNR on 


1 
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resolution and tracking capability, when the spherical signal subspace is used. 
Furthermore, we show that the exponential weighting factor plays a more 
important role than the moving weight. In addition, we see that the resolution 
becomes worse at very low SNR because some noise zeros are approaching to 
the unit circle and can be potentially mistaken for signal zeros. When we use 
the signal subspace without using averaged eigenvalues, we observe that the 
results do not have any significant difference from the case using averaged 
eigenvalues. However, if we use the noise subspace with spherical averaged 
eigenvalues, poor resolution is obtained. Therefore, it is better to adopt the 


signal subspace for resolving two closely-spaced signal sources. 
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TABLE 1. EFFECTIVE STANDARD DEVIATIONS OBTAINED 
BY USING THE TLS ALGORITHM 


SNR (dB) Estimated Deviation for | Estimated Deviation for 
fi=0.23Hz f2=0.27Hz 
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‘ f1=0.23, £2=0.27 
saan > 






30 indépendent trials 
-1.5 7 3 ee ‘ : 
-1.5 -l-0.5 0 0.5 1 iS 
Figure 3.1.1 TLS solution for SNR=20dB, P=8 
N=50, f1=0.23, f2=0.27, undamped sinusoids 


7 : f1=0.23, f2=0.27 


 Phil=pils — 





: 30 independent trials 
-1.5 ; ‘ : ‘ ‘ 
15). t'-05). 05 OS: tt 45 
Figure 3.1.2 TLS solution for SNR=12dB, P=8 
N=50, fi=0.23, f2=0.27, undamped sinusoids 
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: 30 independent trials 





as. 4: SO 10S) ts 
Figure 3.1.3 -TLS solution for SNR=5dB, P=8 
N=50, f1=0.23, f2=0.27, undamped sinusoids 


: £1=0.23, f£2=0.27 
pals phlt=0 


30 independent trials 





5 : : 

-1.5 -1 -0.5 0 0.5 1 1.5 
Figure 3.1.4 TLS solution for SNR=3dB, P=8, 

N=50, f1=0.23, f2=0.27, undamped sinusoids 
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: £1=0.23, f£2=0.27 
aaa es 





30 independent trials 


15 -1 05 O oO5 1 145 
Figure 3.1.5 TLS solution for SNR=1dB, P=8 
N=50, f1=0.23, f2=0.27, undamped sinusoids 








30 independent trials 
“1.5 
13. 1 05. 0 05. °«21~«215 | | 
Figure 3.1.6 TLS solution for SNR=O0dB, P=8 a 


N=50, f1=0.23, f2=0.27, undamped sinusoids 
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: £1=0.25, £2=0.27 





SIS< 2 005°) OOS ke 1S 
Figure 3.2.1 TLS solution for SNR=3dB, P=8 
_N=50, f1=0.25, f2=0.27, undamped sinusoids 


1.5 









: 30 independent trials 






“15 -l 05 O O05 1 15 
Figure 3.2.2 TLS solution for SNR=12dB, P=8 
N=50, f1=0.25, f2=0.27, undamped sinusoids 
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715 -1 0.5 0 05 1 1.5 
Figure 3.2.3 TLS solution for SNR=3dB, P=8. N=50 
f1=0.25, f2=0.27, di=n/4, 62=-1/4, undamped sinusoids 
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“15 -1  -05 0 0.5 1 1.5 
Figure 3.2.4 TLS solution for SNR=12dB, P=8, N=50 
f1=0.25, f2=0.27, o1=1/4, 62=-1/4, undamped sinusoids 
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: £1=0.25, f£2=0:27 
Philepi/4, ae 


: 10 independent trials 


Pg 08) 0 Os, 
Figure 3.3.1 TLS solution for SNR=12dB, P=2 
N=50, f1=0.25, f2=0.27, undamped sinusoids 


: £1=0.25, f2=0.27 
phil=pi/4, ae 


: 10 independent:trials 





“15 -1 O05 O 05 1 15 
Figure 3.3.2 TLS solution for SNR=12dB, P=4 
N=50, f1=0.25, f2=0.27, undamped sinusoids 
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: £1=0.25, f2=0:27 












1.5 a ek ; ; : A 
15 -1 05 O 05 1 15 ‘ 
Figure 3.3.3 TLS solution for SNR=12dB, P=8 a 
N=50, f1=0.25, fa=0.27, undamped sinusoids : 
| 
1.5 . 
'£1=0.25, £2=0.27 | 
-1.5 i . 2 7 ; 
-1.5 -1 -0.5 0 0.5 1 15 
Figure 3.3.4 TLS solution for SNR=12dB, P=18 


N=50, f1=0.25, f£2=0.27, undamped sinusoids 
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TTT $10.25, £2=0.27 
7 a phil =pi4, a 


: 10 independent trials 





Ss a OS Oo OS. 1 2S 
Figure 3.3.5 TLS solution for SNR=12dB, P=28 
N=50, f1=0.25, f2=0.27, undamped sinusoids 


: f1=0.25, £2=0.27 
nto pu 


10 independent trials 





As .t 0S ; 65 1. = 15 
Figure 3.3.6 TLS solution for SNR=12dB, P=38 
N=50, f1=0.25, f2=0.27, undamped sinusoids 


DLs, C 








: f1=0.25, £2=0;27 
pela, areas . 


£0 independent trials 
 destojed signal subspace by: 3 col. 





-1.5 
Figure 3.4.1 TLS solution for SNR=12dB, P=8 a 

N=50, f1=0.25, f2=0.27, undamped sinusoids oe 
signal and noise subspaces mixed by 3 columns : 





#1=0.25, f2=0.27 
phew phi2=0 


Foanieseicn trials 
 desroyed signal subspace by ¢ 6 col. 





mes a O00", 05 OS. OL. AS 
Figure 3.4.2 TLS solution for SNR=12dB, P=8 

N=50, fi=0.25, f2=0.27, undamped sinusoids 
signal and noise subspaces mixed by 6 columns 
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f1=0.25, £2=0.27 


alphal=0.8 ; 
Phlopis phi 


ee 9 


3 30 independent trials 

15 ; : 

-1.5 -1 -0.5 0 0.5 1 15 

Figure 3.5.1 TLS solution for SNR=5dB, P=8 
N=50, f1=0.25, f2=0.27, damped sinusoids 
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alphal=0.8 
aes 9 


: £1=0.25, f£2=0.27 
onan ead 


30 independent trials 
1s. Si - ae 0 0.5 1 1.5 


Figure 3.5.2 TLS solution for SNR=20dB, P=8 
N=50, f1=0.25, f2=0.27, damped sinusoids 
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30 independent trials 


aa: =e 205 ; a5 1 LS 
Figure 3.5.3 TLS solution for SNR=5dB, P=18 
N=50, f1=0.25, f£2=0.27, damped sinusoids 
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: £1=0.25, f2=0.27 
ileal phi2=0 


alphal=0.8 
alpha2=0. 9 


+30 30 independent trials 
15 -1 05 O 0.5 1 1.5 


Figure 3.5.4 TLS solution for SNR=20dB, P=18 
N=50, f1=0.25, f2=0.27, damped sinusoids 
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° yalpha=0. 256 






: oBhjl=pi/4, phi2=0 
pe 3° oving, f2=0. 27 


:10 recursive realizations 


-1.5 -1 05 0 0.5 1 1.5 
Figure 4.1.1 RTLS solution for SNR=12dB,P=8,N=100 
a=0.25, weight=3%, f1 is time-varying, using 
the signal subspace information & performing 
the eigendecomposition of spherical matrix 'S' 


| 
3 
: 
| 





0 2 4 6 8 10 
Figure 4.1.2 averaged eigenvalues of updated spherical matrix 'S' 
SNR=12dB,P=8,N=100,a0=0.25,weight=3%,f1 is time-varying 
solution derived by using the signal subspace information 
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‘alpha=0.95 : phil=pi/4, phi2=0 
: weight=0.03 : 


: 10 independent trials : 
: 10 recursive realizations 
15 -1 05 #0 O5 1. 15 
Figure 4.1.3 RTLS solution for SNR=12dB,P=8,N=100 
a=0.95, weight=3%, fi is time-varying, using 
the signal subspace information & performing 
the eigendecomposition of spherical matrix 'S' 
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300 





0 2 4 6 8 10 
Figure 4.1.4 averaged eigenvalues of updated spherical matrix 'S' 
SNR=12dB,P=8,N=100,a=0.95 ,weight=3%,f1 is time-varying 
solution derived by using the signal subspace information 
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phil =pi/4, phi2=0 
fl moving, f2=0. 27 






: alpha=0.95 
ste 10: 











40 independent trials * 
: 10 recursive realizations 


1s i as OOS aS 
Figure 4.2.1 RTLS solution for SNR=12dB,P=8,N=100 
a=0.95, weight=10%, f1 is time-varying, using 
the signal subspace information & performing 
the eigendecomposition of spherical matrix 'S' 





° 2 4 6 8 10 
Figure 4.2.2 averaged eigenvalues of updated spherical matrix 'S' 
SNR=12dB,P=8,N=100,0=0.95,weight=10%,f1 is time-varying 
solution derived by using the signal subspace information 
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:alpha=0.95 : phil=pi/4, phi2=0 
hace ae 10 : fl a f2=0. 27 






‘YO independeni trials : 
‘10 recursive realizations 
-1.5 -1 -0.5 0 0.5 1 1.5 
Figure 4.3.1 RTLS solution for SNR=12dB, P=5, N=100 © 
a=0.95, weight=10%, f1 is time-varying, using 
the signal subspace information & performing 
the eigendecomposition of spherical matrix 'S' 
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: alpha=0.95 


: phil=pi/4, phi2=0 
: weight=0. 10 : 


fl Lmoyne: P=0. 27 






10 Independent trials : 
:10 recursive realizations 
-1.5 -1 ry 3 0 0.5 1 LS 
Figure 4.3.2 RTLS solution for SNR=12dB, P=3, N=100 
a=0.95, weight=10%, f1 is time-varying, using 
the signal subspace information & performing 
the eigendecomposition of spherical matrix 'S' 
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: alpha=0. 95 
Nae 10 : 


: phil=pi/4, phi2=0 
f1 moving, £2=0. 21 






10 ok trials : 


:10 recursive realizations 


-1.5 : 
-1.5 -1 -0.5 0 0.5 1 1.5 


Figure 4.4.1 RTLS solution for SNR=0dB,P=8,N=100 
a=0.95, weight=10%, f1 is time-varying, using 
the signal subspace information & performing 
the eigendecomposition of spherical matrix 'S' 
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:alpha=0.95  : phil=pi/4, phi2=0 
 welen =0,10 sg esie ee f2=0.27 





:10 independent trials : 
:10 recursive realizations 


“15 ll 05.0. 05. +1 415 
Figure 4.4.2 RTLS solution for SNR=-5dB,P=8,N=100 
a=0.95, weight=10%, f1 is time-varying, using 
the signal subspace information & performing 
the eigendecomposition of spherical matrix 'S' 
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:alpha=0.95  ; phil=pi/4, phi2=0 
a Ts 10 : fl mys f2=0. 27 


10 independent trials : 
10 recursive realizations 


5. 1°05 0 05. +1 415 
Figure 4.5.1 RTLS solution for SNR=12dB,P=8,N=100 
a=0.95, weight=10%, fi is time-varying, using 
the signal subspace information & performing 
the eigendecomposition of spherical matrix 'R' 
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0 Z 4 6 8 10 
Figure 4.5.2 non-averaged eigenvalues of correlation matrix 'R' 
SNR=12dB,P=8,N=100,a=0.95,weight=10%,f1 is time-varying 
solution derived by using the signal subspace information 
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:alpha=0.95 
: weight=0.10 : 





-1.5 -1 -0.5 0 0.5 1 1.5 
Figure 4.6.1 RTLS solution for SNR=0dB,P=8,N=100 
a=0.95, weight=10%, f1 is time-varying, using 
the signal subspace information & performing 
the eigendecomposition of spherical matrix 'R' 
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0 2 4 6 8 10 
Figure 4.6.2 non-averaged eigenvalues of correlation matrix 'R' 
SNR=0dB,P=8,N=100,0.=0.95, weight=10%,f1 is time-varying 
solution derived by using the signal subspace information 
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Mo independent trials 
:10 recursive realizations 


3. a 05 0 05 1 ~«215 
Figure 4.6.3. RTLS solution for SNR=-5dB ,P=8,N=100 
a=0.95, weight=10%, fi is time-varying, using 
the signal subspace information & performing 
the eigendecomposition of spherical matrix 'R' 
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0 Z 4 6 8 10 
Figure 4.6.4 non-averaged eigenvalues of correlation matrix 'R’ re oe 
SNR=-5dB,P=8,N=100,a=0.95,weight=10%,f1 is time-varying 
solution derived by using the signal subspace information 
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"10 independent trials 
:10 recursive realizations 


-1.5 : 
-1.5 -1 “0. 5 0 0.5 1 1.5 


Figure 4.7.1 RTLS solution for SNR=12dB,P=8,N=100 
=0.95, weight=10%, fi is time-varying, using 
the noise subspace information & performing 
the eigendecomposition of spherical matrix 'S' 


0 2 4 6 8 10 
Figure 4.7.2 averaged eigenvalues of updated spherical matrix 'S' 
SNR=12dB,P=8,N=100,a=0.95,weight=10%,f1 is time-varying 
solution derived by using the noise subspace information 
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10 Greig geiacat trials 
: 10 recursive realizations 
rs 
15-1 05. 0 O05 1 45 , ‘ 
Figure 4.8.1 RTLS solution for SNR=12dB,P=5,N=100 — ‘3 
a=0.95, weight=10%, f1 is time-varying, using . mu 
the noise subspace information & performing 
the eigendecomposition of spherical matrix 'S' 
1.5 
:alpha=0.95  : phil=pi/4, phi2=0 
acts 10 : fil moving, f2=0. 27 
1 cere erage hae aE Or errr 
yy Rene eee ae corr erties mene opener Cosy es 
Ol eske ie Ana ern at ea gee 
OVS es cere at sh acs aac ana oan xan 
SP Lite ee a kee ere 
10 independent tials : 
:10 recursive realizations 
-1.5 . . . . . 
-1.5 -1 -0.5 0 0.5 1 1.5 
Figure 4.8.2 RTLS solution for SNR=12dB ,P=3,N=100 fe 
a=0.95, weight=10%, f1 is time-varying, using 
the noise subspace information & performing 


the eigendecomposition of spherical matrix 'S' 
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id independent trials : 
:10 recursive realizations 
= aS 0. 5 0 05 1 15 
Figure 4.9.1 RTLS solution for SNR=15dB ,P=8,N=100 
a=0.95, weight=10%, fi is time-varying, using 
the noise subspace information & performing 
the eigendecomposition of spherical matrix 'S' 
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calpha=9.95 phil=pi/4, phi2=0 
weight. 10: fl madving, £2=0. va 









40 Madehendieat ‘rials 
10 recursive realizations 


21.5 -1 0. 5 0 0.5 1 iS 
Figure 4.9.2 RTLS solution for SNR=30dB,P=8,N=100 
a=0.95, weight=10%, f1 is time-varying, using 
the noise subspace information & performing 
the eigendecomposition of spherical matrix 'S' 
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APPENDIX A : MATLAB SOURCE CODE 


This appendix consists of MATLAB codes used to implement the Total 
Least Squares ( TLS ) algorithm mentioned in Chapter III and the Recursive 


Total Least Squares ( RTLS ) algorithm mentioned in Chapter IV. 


1. IMPLEMENTATION OF THE TLS ALGORITHM 
clear 

Yinitialize the set of roots 

rt=[]; 

nn=input(' trial number '); 

SNR=input(' S/N(dB) = '); 

p=input(' prediction filter order '); 


y 1=zeros(p,1);y2=zeros(p,1); 





Joconstruct the received signal and begin the trial 
for h=1:nn, 

n=0:49; 

rand(‘normal'); 


wl=rand(1,50); 


rand(‘normal’); 

w2=rand(1,50); 

f1=0.25; %signal frequencies 
f2=0.27; %f1=0.23 for option 
phil=pi/4; %signal phase angles 
phi2=0; | 
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alphal=0.8; 
alpha2=0.9; 
j=sqrt(-1); 
wh=w 1+j*w2; 
sl=exp(j*(2*pi*fl *n+phil )); 
s2=exp(j*(2*pi*f2*n+phi2)); 
sl=exp(-alphal] *n+j*(2*pi*fl *n+phil)); 
s2=exp(-alpha2*n+j*(2*pi*f2*n+phi2)), 
ss=sl+s2; 
xn=ss+sqrt(1/(2*10(SNR/10))).*wn; 
nsin=2; 
n=50; 
if n>p+2, 

for k=] :n-p, 

A(k, 1:p)=xn(k:p+k-1); 
end 
b=xn(1,p+1:n); 
bb=conj(b)’; 

end . 
C=[bb A]; 
[u,s,vJ=svd(C); 
sigma=diag(s); 
ct=v(1,nsin+1:p+1); 
c=conj(ct)’; 


vp=v(2:p+1 ,nsin+1:p+1); 
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%decaying constants 


%independent complex noise 


%create two undamped signals 


%create two damped signals 


%received resulted signal 


%for two sinusoidal sources 


%window data sample number 


JYogenerate data matrix 


Jogenerate observation vector 


Yocreate augmented matrix 


Yoperform singular value decom. 


%extract the desired vector 


v2=[ctivp]; 
xtls=-(vp*conj(c))./(conj(ct)*c); 
z=[1 -fliplr(conj(xtls)')]; 
rt=[rt roots(z)]; 
z1= 0.3+j*0.7; 
zZ2=-0.3+j*0.7; 
real_rt=real(rt); 
for nl=1:p, 
if real_rt(n1,h) > 0, 
yl(nl1,h)=rt(n1,h); 
yll1=sort(y1); 
yl1l=flipud(y11); 


new_y11(1:p/2,h)=y11(1:p/2,h); 





%solution from noise subspace 


Yoprediction filter polynomial 


Yocreate sector regions vs. origin 


%calculate estimated deviations 


k1=find((angle(new_y11(1:p/2,h)) > angle(z1))... 
& (angle(new_y11(1:p/2,h)) < pi/2)); 


if length(k1)~=1, 


error('More than 1 element of new_yl1 is nonzero.’), 


end 


esti_y11(1,h)=new_y11(k1,h); 


estd1=sqrt(1/nn*(sum(abs(esti_y11(1,h)... 
-(0.1253+j*0.9921) .*ones(1,nn)) .42))); 


else 
y2(n1,h)=rt(n1 ,h); 
y22=sort(y2); 
y22=flipud(y22); 





new_y22(1:p/2,h)=y22(1:p/2,h); 
k2=find((angle(new_y22(1:p/2,h)) > pi/2)... 
& (angle(new_y22(1:p/2,h)) < angle(z2))); 
if length(k2)~=1, 
error(More than 1 element of new_y22 is nonzero.'); 
end 
esti_y22(1,h)=new_y22(k2,h); 
estd2=sqrt(1/nn*(sum(abs(esti_y22(1,h)... 
-(-0.1253+j*0.9921) .*ones(1,nn)) .42))); 
end 
end 

end 
%plot the zeros distribution situations 
axis(‘square’) 
axis([-1.5 1.5 -1.5 1.5]) 
plot(sin(0:2*pi/100:2*pi),cos(0:2*pi/100:2*pi),'-w’,... 
real(new_y11),imag(new_y11),'o',real(new_y22),imag(new_y22),'o'); 
grid | 
title TLS solution for SNR= , P= , N= '); 
text(-1.1,1.3,' alphal= '); 
text(-1.1,1.1,' alpha2= '); 
text(0,1.3,' fl= , f2= '); 
text(0,1.1,’ phil= , phi2= '); 
text(0,-1.3,' 30 independent trials '); 


text(-1,-1.4,' destroyed signal subspace by # columns ' ); 
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JREMARKS: 

%(1).When SNR down to 1 dB, we need change z2 to -0.4+j*0.8 
% in order to make k2 for new_y22 satisfy ‘find’ condition. 
%(2).For SNR down to 0 dB, we have to change z]1 to 0.4+j*0.7 


% in order to make k1 for new_y11 satisfy ‘find’ condition. 
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2. IMPLEMENTATION OF THE RTLS ALGORITHM USING 
THE SIGNAL SUBSPACE AND PERFORMING THE EIGEN- 
DECOMPOSITION OF MATRIX ' S' 

Note that S = a*D+(1-c)*B*B# B = Ulex 

clg 

clear 

%initialization of conditions 

tt=[];esp=[];sig=[];noi=[]; 

nt=input(' trial number= '); 

ni=input(’ iteration number= '); 

p =input(' prediction order= '); 

SNR=input(' signal-to-noise ratio, S/N(dB)= '); 

weight=input(' the weight of moving frequency= '); 

alpha=input(' exponential weighting factor, alpha= '); 

Y%construct received signal and begin independent trial 

for h=1:nt, 

nl=1:100; 

rand(‘normal'); 

wl=rand(1,100); 

rand(‘normal’); 

w2=rand(1,100); 

f1=0.23-(h-1)*weight/nt; . Yomake fi time-varying moving 
£2=0.27; | 

phil=pi/4; 

phi2=0; 


TL, ots 


j=sqrt(-1); 
wn=W1+j*w2; 
sl=exp(j*(2*pi*fl *n1+phi1)); 
s2=exp(j*(2*pi*f2*n1+phi2)); 
ss=sl+s2; . 
xn=ss+sqrt(1/(2*104(SNR/10))).*wn; 
nsin=2; 
n2=50; 
if n2>p+2, 
for k=1:n2-p, 
A(k, 1 :p)=xn(k:p+k- 1); 
end 
b=xn(1,p+1:n2); 
bb=conj(b)’; 
end 
B=[bb A]; 
aph=alpha .4n1; 
aph=flipIr(aph(1:n2-p)); 
T=sqrt(diag(aph)); 
C=FtB; 
R=C™*C; 
[U,D]=eig(R); 
for t=I:ni, 
y=xn(1,n2-p+t:n2+t); 


—y" 
X=y, 


te 





Yindependent complex noise 


%create two undamped signals 


Jreceived resulted signal 
%for two sinusoidal sources 


Jowindow data sample number 


%generate data matrix 


JYgenerate observation vector 


JYocreate augmented matrix 
%define estd. correlation matrix 
%perform eigendecomposition 


%recursively update subspace 








beta=U'*x; 

D1=sort((diag(D))); 
D2=flipud(D1); 

D2=D2'; 

s_eigval=mean(D2(1 :nsin)); 

n_eigval=mean(D2(nsin+1:p+1)); 

s_eigspa=s_eigval .*eye(nsin); 

n_eigspa=n_eigval .*eye(p-nsin+1); 

s_diaspa=diag(s_eigspa); 

n_diaspa=diag(n_eigspa); 

Dhat1l=[s_eigspa zeros(nsin,p-nsin+1)); 
Dhat2=[zeros(p-nsin+1,nsin) n_eigspa]; 

Dhat =[Dhat1;Dhat2]; 

S=alpha .*Dhat+(1-alpha) .*beta*beta’;, %our major concern 
[UU,DD)=eig(S); 

Unew=U*UU; 

Dnew=DD; 

gt=Unew(1,l:nsin); | | ; %extract the desired vector 
g=conj(gt)’; 

ulp=Unew(2:p+1,1:nsin); 

ul=[gt;ulp]; 

xtls=(ulp*conj(g))./(1-conj(gt)*g); %solution from signal subspace 
z=[1 -fliplr(conj(xtls)’)]; 

rt=[rt roots(z)]; 


esp=[esp D2']; 








sig=[sig s_diaspa]; 
noi=[noi n_diaspa)J; 
U=Unew; Jassign new U & D for updating 
D=Dnew; 
end 
end 
%plot zeros in complex z-plane 
axis(‘square'); 
axis({-1.5 1.5 -1.5 1.5]); 
plot(sin(0:2*pi/100:2*pi),cos(0:2*pi/100:2*pi),'-w’,... 
real(rt),imag(rt),'0'), grid 
title’ RTLS solution for SNR= dB, P= , N= '); | 
text(-1,1.3,' alpha= ‘); 
text(-1,1.1,’ weight= '); 
text(0.1,1.1,' {1 moving, f2=0.27 '); 
text(0.1,1.3,' phil=pi/4, phi2=0'); 
text(-0.5,-1.2,' 10 independent trials °); 
text(-0.5,-1.4,' 10 recursive realizations '); 
%plot spherical singularvalues status of updated matrix 
clg 
axis; 
n3=1:p+1; 
plot(n3,esp,'o'),grid 
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- 3. IMPLEMENTATION OF THE RTLS ALGORITHM USING 
THE SIGNAL SUBSPACE AND PERFORMING THE EIGEN- 
DECOMPOSITION OF MATRIX 'R' 

Note that R = a*R+(1-c)*x*xH 

clear 
Joinitialization of eonditions 
rt=[];esp=[]; 
nt=input(’ trial number= '); 
ni=input(’ iteration number= '); 
p =input(' prediction order= ‘); 
SNR=input(' signal-to-noise ratio, S/N(dB)= '); 
weight=input(' the weight of moving frequency= '); 
alpha=input(' exponential weighting factor, alpha= ‘); 
Joconstruct received signal and begin independent trial 
for h=1:nt, 

nl=1:100; 

rand(‘normal'); 

w 1=rand(1,100); 

rand(‘normal'); 

w2=rand(1,100); 

£1=0.23-(h-1)*weight/nt; | %make fi time-varying moving 

f2=0.27; 

phil=pi/4; 

phi2=0; 

j=sqrt(-1); 
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wn=w1+j*w2; 
sl=exp(j*(2*pi*f1*n1+phil)); 
s2=exp(j*(2*pi*f2*n1+phi2)); 
ss=sl+s2; 
xn=ss-+sqrt(1/(2*104(SNR/10))).*wn; 
nsin=2; 
n2=50; 
if n2>p+2, 
for k=1:n2-p, 
A(k,1:p)=xn(k:p+k-1); 
end 
=xn(1,p+1:n2); 
bb=conj(b)’; 
end 
B=[bb A]; 
aph=alpha .4n1; 
aph=fliplr(aph(1 :n2-p)); 
T=sqrt(diag(aph)); 
C=T*B; 
R=C'*C; 
[U,D]=eig(R),; 
for t=1:ni, 
y=xn(1,n2-p+t:n2+t); 
x=y'; 


R=alpha .*R+(1-alpha) y's 
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%independent complex noise 


Yocreate two undamped signals 


%received resulted signal 
%for two sinusoidal sources 


Yowindow data sample number 


Y%generate data matrix ' 


J%generate observation vector 


Yocreate augmented matrix 
%define estd. correlation matrix 
%perform eigendecomposition 


Yrecursively update subspace 


Your major concern 








[UU,DD]=eig(R); 
Di=flipud(sort((diag(DD)))); _ | 
gt=UU(1,1 :nsin); %extract the desired vector 
g=conj(gt)'; | 
ulp=UU(2:p+1,1:nsin); 
ul=[gt;ulp]; 
xtls=(ul p*conj(g))./(1-conj(gt)*g); ~~ %solution from signal subspace 
z=([1 -fliplr(conj(xtls)')]; 
rt=[rt roots(z)]; 
esp=[esp D1]; 
end 
end 
Yplot zeros in complex z-plane 
axis(‘square’); | 
axis({-1.5 1.5 -1.5 1.5]); 
plot(sin(0:2*pi/100:2*pi),cos(0:2*pi/100:2*pi),'-w’,... 
real(rt),imag(rt),'o'),gtid 
title RTLS solution for SNR= dB, P= , N= '); 
text(-1,1.3,' alpha= ‘); 
text(-1,1.1,' weight= '); 
text(0.1,1.1,' f1 moving, f2=0.27 '); 
text(0.1,1.3,' phil=pi/4, phi2=0'); 
text(-0.5,-1.2,' 10 independent trials "); 
text(-0.5,-1.4,' 10 recursive realizations '); 


clg 


a 


axis; 
n3=1:p+1; 
plot(n3,esp,'o’),grid 
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4. IMPLEMENTATION OF THE RTLS ALGORITHM USING 
THE NOISE SUBSPACE AND PERFORMING THE EIGEN- 
DECOMPOSITION OF MATRIX ' S ' | 
Note that S$ = a*D+(1-)*B*Bo B = UH+x 

clg 

clear 

%initialization of conditions 

rt=[];esp=[];sig=[J;noi=[]; 

nt=input(' trial number= '); 

ni=input(' iteration number= '); 

p =input(' prediction order= '); 

SNR=input(' signal-to-noise ratio, S/N(dB)= '); 

weight=input(’ the weight of moving frequency= '); 

alpha=input(' exponential weighting factor, alpha= '); 

Yconstruct received signal and begin independent trial 

for h=1:nt, 

~ n1=1:100; 
rand(‘normal'’); 
wl=rand(1,100); 
rand(‘normal'’); 
w2=rand(1,100); 

f1=0.23-(h-1)*weight/nt; make fi time-varying moving 
£2=0.27: 

phil=pi/4; 

phi2=0; 
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j=sqrt(-1); 
wn=w1+j*w2; 
sl=exp(j*(2*pi*f1*nl+phil)); 
s2=exp(j*(2*pi*f2*n1+phi2)); 
ss=sl+s2; 
xn=ss+sqrt(1/(2*104(SNR/10))).*wn; 
nsin=2; 
n2=50; 
if n2>p+2, 
for k=1:n2-p, 
A(k,1:p)=xn(k:p+k- 1); 
end 
b=xn(1,p+1:n2); 
bb=conj(b)’; 
end 
B=[bb AJ; 
aph=alpha .4n1; 
aph=fliplr(aph(1:n2-p)); 
T=sqrt(diag(aph)), 
C=T*B; 
R=C'*C; 
[U,D]=eig(R); 
for t=1:ni, 
y=xn(1,n2-p+t:n2+t); 


' 
xX=y; 
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J%independent complex noise 


Yocreate two undamped signals 


%received resulted signal 
%for two sinusoidal sources 


%window data sample number 


Yogenerate data matrix 


%generate observation vector 


%create augmented matrix 
%define estd. correlation matrix 
%perform eigendecomposition 


%recursively update subspace 





beta=U'*x; 
D1=sort((diag(D))); 
D2=flipud(D1); 
D2=D2'; 
s_eigval=mean(D2(1 :nsin)); 
n_eigval=mean(D2(nsin+1:p+1)); 
s_eigspa=s_eigval .*eye(nsin); 
n_eigspa=n_eigval .*eye(p-nsin+1); 
s_diaspa=diag(s_eigspa), 
n_diaspa=diag(n_eigspa); 
Dhati=[s_eigspa zeros(nsin,p-nsin+1)]; 
Dhat2=[zeros(p-nsin+1,nsin) n_eigspa]; 
Dhat =[Dhat1;Dhat2]}; 
S=alpha .*Dhat+(1-alpha) .*beta*beta’; %our major concern 
[UU,DD]=eig(S); 
Unew=U*UU; 
Dnew=DD; 
ct=Unew(1,nsin+1:p+1); %extract the desired vector 
c=conj(ct)’; 
vp=Unew(2:p+1,nsin+1:p+1); 
v2=[ct;vp]; . | 
xtls=-(vp*conj(c))./(1-conj(ct)*c); %solution from noise subspace 
z=[1 -flipir(conj(xts)’)]; 
tt=[rt roots(z)]; 


esp=[esp D2']; 
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| sig=[sig s_diaspa]; 
noi=(noi n_diaspa]; 
U=Unew; %assign new U & D for updating 
D=Dnew; 
end 
end 
J%plot zeros in complex z-plane 
axis('square’); 
axis([-1.5 1.5 -1.5 1.5]); 
plot(sin(0:2*pi/100:2*pi),cos(0:2*pi/100:2*pi),'-w',... 
real(rt),imag(rt),'o'),grid 
title’ RTLS solution for SNR= dB, P= , N= '); 
text(-1,1.3,’ alpha= '); 
text(-1,1.1,' weight= '); 
text(0.1,1.1,' f1 moving, f2=0.27 '); 
text(0.1,1.3,’ phil=pi/4, phi2=0'); 
text(-0.5,-1.2,’ 10 independent trials ‘); 
text(-0.5,-1.4,' 10 recursive realizations '); 
Yplot spherical singularvalues status of updated matrix 
clg 
axis; 
n3=1:p+1; 
plot(n3,esp,'o'),grid 
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